home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / curve_isect / Bezier.cc next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  12.4 KB  |  368 lines
  1. #include <stddef.h> // for size_t
  2. // #include <sys/stdtypes.h> // for size_t, on some systems
  3. #include "Bezier.h"
  4. #include <math.h>
  5. /* The value of 1.0 / (1L<<23) is float machine epsilon. */
  6. #ifdef FLOAT_ACCURACY
  7. #define INV_EPS (1L<<23)
  8. #else
  9. /* The value of 1.0 / (1L<<14) is enough for most applications */
  10. #define INV_EPS (1L<<14)
  11. #endif
  12. #define log2(x) (log(x)/log(2.))
  13.  
  14. extern "C" void qsort( char *base, int nel, size_t width, int (*compar)(void *, void *));
  15.  
  16. int compare_doubles( void *a, void *b )
  17.     {
  18.     register double *A = (double *)a, *B = (double *)b;
  19.     return( *A > *B )?1:(*A < *B ? -1 : 0 );
  20.     }
  21.  
  22. void Sort( double *array, int length )
  23.     {
  24.     qsort( (char *)array, length, sizeof( double ), compare_doubles );
  25.     }
  26.     
  27. /*
  28.  * Split the curve at the midpoint, returning an array with the two parts
  29.  * Temporary storage is minimized by using part of the storage for the result
  30.  * to hold an intermediate value until it is no longer needed.
  31.  */
  32. #define left r[0]
  33. #define right r[1]
  34. Bezier *Bezier::Split()
  35.     {
  36.     Bezier *r = new Bezier[2];
  37.     (left.p0 = p0)->refcount++;
  38.     (right.p3 = p3)->refcount++;
  39.     left.p1 = new point(  p0, p1, 1);
  40.     right.p2 = new point( p2, p3, 1);
  41.     right.p1 = new point( p1, p2, 1); // temporary holding spot
  42.     left.p2 = new point ( left.p1, right.p1, 1);
  43.     *right.p1 = mid( right.p1, right.p2 ); // Real value this time
  44.     left.p3 = right.p0 = new point( left.p2, right.p1, 2 );
  45.     return r;
  46.     }
  47. #undef left
  48. #undef right
  49.  
  50.     
  51. /*
  52. * Test the bounding boxes of two Bezier curves for interference.
  53. * Several observations:
  54. *    First, it is cheaper to compute the bounding box of the second curve
  55. *    and test its bounding box for interference than to use a more direct
  56. *    approach of comparing all control points of the second curve with 
  57. *    the various edges of the bounding box of the first curve to test
  58. *     for interference.
  59. *    Second, after a few subdivisions it is highly probable that two corners
  60. *    of the bounding box of a given Bezier curve are the first and last 
  61. *    control point.  Once this happens once, it happens for all subsequent
  62. *    subcurves.  It might be worth putting in a test and then short-circuit
  63. *    code for further subdivision levels.
  64. *    Third, in the final comparison (the interference test) the comparisons
  65. *    should both permit equality.  We want to find intersections even if they
  66. *    occur at the ends of segments.
  67. *    Finally, there are tighter bounding boxes that can be derived. It isn't
  68. *    clear whether the higher probability of rejection (and hence fewer
  69. *    subdivisions and tests) is worth the extra work.
  70. */
  71. int IntersectBB( Bezier a, Bezier b )
  72.     {
  73.     // Compute bounding box for a
  74.     double minax, maxax, minay, maxay;
  75.     if( a.p0->x > a.p3->x )     // These are the most likely to be extremal
  76.     minax = a.p3->x, maxax = a.p0->x;
  77.     else
  78.     minax = a.p0->x, maxax = a.p3->x;
  79.     if( a.p2->x < minax )
  80.     minax = a.p2->x;
  81.     else if( a.p2->x > maxax )
  82.     maxax = a.p2->x;
  83.     if( a.p1->x < minax )
  84.     minax = a.p1->x;
  85.     else if( a.p1->x > maxax )
  86.     maxax = a.p1->x;
  87.     if( a.p0->y > a.p3->y )         
  88.     minay = a.p3->y, maxay = a.p0->y;
  89.     else
  90.     minay = a.p0->y, maxay = a.p3->y;
  91.     if( a.p2->y < minay )
  92.     minay = a.p2->y;
  93.     else if( a.p2->y > maxay )
  94.     maxay = a.p2->y;
  95.     if( a.p1->y < minay )
  96.     minay = a.p1->y;
  97.     else if( a.p1->y > maxay )
  98.     maxay = a.p1->y;
  99.     // Compute bounding box for b
  100.     double minbx, maxbx, minby, maxby;
  101.     if( b.p0->x > b.p3->x )         
  102.     minbx = b.p3->x, maxbx = b.p0->x;
  103.     else
  104.     minbx = b.p0->x, maxbx = b.p3->x;
  105.     if( b.p2->x < minbx )
  106.     minbx = b.p2->x;
  107.     else if( b.p2->x > maxbx )
  108.     maxbx = b.p2->x;
  109.     if( b.p1->x < minbx )
  110.     minbx = b.p1->x;
  111.     else if( b.p1->x > maxbx )
  112.     maxbx = b.p1->x;
  113.     if( b.p0->y > b.p3->y )         
  114.     minby = b.p3->y, maxby = b.p0->y;
  115.     else
  116.     minby = b.p0->y, maxby = b.p3->y;
  117.     if( b.p2->y < minby )
  118.     minby = b.p2->y;
  119.     else if( b.p2->y > maxby )
  120.     maxby = b.p2->y;
  121.     if( b.p1->y < minby )
  122.     minby = b.p1->y;
  123.     else if( b.p1->y > maxby )
  124.     maxby = b.p1->y;
  125.     // Test bounding box of b against bounding box of a
  126.     if( ( minax > maxbx ) || ( minay > maxby )  // Not >= : need boundary case
  127.     || ( minbx > maxax ) || ( minby > maxay ) )
  128.     return 0; // they don't intersect
  129.     else
  130.     return 1; // they intersect
  131.     }
  132.     
  133. /* 
  134. * Recursively intersect two curves keeping track of their real parameters 
  135. * and depths of intersection.
  136. * The results are returned in a 2-D array of doubles indicating the parameters
  137. * for which intersections are found.  The parameters are in the order the
  138. * intersections were found, which is probably not in sorted order.
  139. * When an intersection is found, the parameter value for each of the two 
  140. * is stored in the index elements array, and the index is incremented.
  141. * 
  142. * If either of the curves has subdivisions left before it is straight
  143. *    (depth > 0)
  144. * that curve (possibly both) is (are) subdivided at its (their) midpoint(s).
  145. * the depth(s) is (are) decremented, and the parameter value(s) corresponding
  146. * to the midpoints(s) is (are) computed.
  147. * Then each of the subcurves of one curve is intersected with each of the 
  148. * subcurves of the other curve, first by testing the bounding boxes for
  149. * interference.  If there is any bounding box interference, the corresponding
  150. * subcurves are recursively intersected.
  151. * 
  152. * If neither curve has subdivisions left, the line segments from the first
  153. * to last control point of each segment are intersected.  (Actually the 
  154. * only the parameter value corresponding to the intersection point is found).
  155. *
  156. * The apriori flatness test is probably more efficient than testing at each
  157. * level of recursion, although a test after three or four levels would
  158. * probably be worthwhile, since many curves become flat faster than their 
  159. * asymptotic rate for the first few levels of recursion.
  160. *
  161. * The bounding box test fails much more frequently than it succeeds, providing
  162. * substantial pruning of the search space.
  163. *
  164. * Each (sub)curve is subdivided only once, hence it is not possible that for 
  165. * one final line intersection test the subdivision was at one level, while
  166. * for another final line intersection test the subdivision (of the same curve)
  167. * was at another.  Since the line segments share endpoints, the intersection
  168. * is robust: a near-tangential intersection will yield zero or two
  169. * intersections.
  170. */
  171. void RecursivelyIntersect( Bezier a, double t0, double t1, int deptha,
  172.                Bezier b, double u0, double u1, int depthb,
  173.                double **parameters, int &index )
  174.     {
  175.     if( deptha > 0 )
  176.     {
  177.     Bezier *A = a.Split();
  178.     double tmid = (t0+t1)*0.5;
  179.     deptha--;
  180.     if( depthb > 0 )
  181.         {
  182.         Bezier *B = b.Split();
  183.         double umid = (u0+u1)*0.5;
  184.         depthb--;
  185.         if( IntersectBB( A[0], B[0] ) )
  186.         RecursivelyIntersect( A[0], t0, tmid, deptha,
  187.                       B[0], u0, umid, depthb,
  188.                       parameters, index );
  189.         if( IntersectBB( A[1], B[0] ) )
  190.         RecursivelyIntersect( A[1], tmid, t1, deptha,
  191.                       B[0], u0, umid, depthb,
  192.                       parameters, index );
  193.         if( IntersectBB( A[0], B[1] ) )
  194.         RecursivelyIntersect( A[0], t0, tmid, deptha,
  195.                       B[1], umid, u1, depthb,
  196.                       parameters, index );
  197.         if( IntersectBB( A[1], B[1] ) )
  198.         RecursivelyIntersect( A[1], tmid, t1, deptha,
  199.                       B[1], umid, u1, depthb,
  200.                       parameters, index );
  201.         }
  202.     else
  203.         {
  204.         if( IntersectBB( A[0], b ) )
  205.         RecursivelyIntersect( A[0], t0, tmid, deptha,
  206.                       b, u0, u1, depthb,
  207.                       parameters, index );
  208.         if( IntersectBB( A[1], b ) )
  209.         RecursivelyIntersect( A[1], tmid, t1, deptha,
  210.                       b, u0, u1, depthb,
  211.                       parameters, index );
  212.         }
  213.     }
  214.     else
  215.     if( depthb > 0 )
  216.         {
  217.         Bezier *B = b.Split();
  218.         double umid = (u0 + u1)*0.5;
  219.         depthb--;
  220.         if( IntersectBB( a, B[0] ) )
  221.         RecursivelyIntersect( a, t0, t1, deptha,
  222.                       B[0], u0, umid, depthb,
  223.                       parameters, index );
  224.         if( IntersectBB( a, B[1] ) )
  225.         RecursivelyIntersect( a, t0, t1, deptha,
  226.                       B[0], umid, u1, depthb,
  227.                       parameters, index );
  228.         }
  229.     else // Both segments are fully subdivided; now do line segments
  230.         {
  231.         double xlk = a.p3->x - a.p0->x;
  232.         double ylk = a.p3->y - a.p0->y;
  233.         double xnm = b.p3->x - b.p0->x;
  234.         double ynm = b.p3->y - b.p0->y;
  235.         double xmk = b.p0->x - a.p0->x;
  236.         double ymk = b.p0->y - a.p0->y;
  237.         double det = xnm * ylk - ynm * xlk;
  238.         if( 1.0 + det == 1.0 )
  239.         return;
  240.         else
  241.         {
  242.         double detinv = 1.0 / det;
  243.         double s = ( xnm * ymk - ynm *xmk ) * detinv;
  244.         double t = ( xlk * ymk - ylk * xmk ) * detinv;
  245.         if( ( s < 0.0 ) || ( s > 1.0 ) || ( t < 0.0 ) || ( t > 1.0 ) )
  246.             return;
  247.         parameters[0][index] = t0 + s * ( t1 - t0 );
  248.         parameters[1][index] = u0 + t * ( u1 - u0 );
  249.         index++;
  250.         }
  251.         }
  252.     }
  253.  
  254. inline double log4( double x ) { return 0.5 * log2( x ); }
  255.     
  256. /*
  257.  * Wang's theorem is used to estimate the level of subdivision required,
  258.  * but only if the bounding boxes interfere at the top level.
  259.  * Assuming there is a possible intersection, RecursivelyIntersect is
  260.  * used to find all the parameters corresponding to intersection points.
  261.  * these are then sorted and returned in an array.
  262.  */
  263.  
  264. double ** FindIntersections( Bezier a, Bezier b )
  265.     {
  266.     double **parameters = new double *[2];
  267.     parameters[0] = new double[9];
  268.     parameters[1] = new double[9];
  269.     int index = 0;
  270.     if( IntersectBB( a, b ) )
  271.     {
  272.     vector la1 = vabs( ( *(a.p2) - *(a.p1) ) - ( *(a.p1) - *(a.p0) ) );
  273.     vector la2 = vabs( ( *(a.p3) - *(a.p2) ) - ( *(a.p2) - *(a.p1) ) );
  274.     vector la;
  275.     if( la1.x > la2.x ) la.x = la1.x; else la.x = la2.x;
  276.     if( la1.y > la2.y ) la.y = la1.y; else la.y = la2.y;
  277.     vector lb1 = vabs( ( *(b.p2) - *(b.p1) ) - ( *(b.p1) - *(b.p0) ) );
  278.     vector lb2 = vabs( ( *(b.p3) - *(b.p2) ) - ( *(b.p2) - *(b.p1) ) );
  279.     vector lb;
  280.     if( lb1.x > lb2.x ) lb.x = lb1.x; else lb.x = lb2.x;
  281.     if( lb1.y > lb2.y ) lb.y = lb1.y; else lb.y = lb2.y;
  282.     double l0;
  283.     if( la.x > la.y ) 
  284.         l0 = la.x;
  285.     else 
  286.         l0 = la.y;
  287.     int ra;
  288.     if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 ) 
  289.         ra = 0;
  290.     else
  291.         ra = (int)ceil( log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
  292.     if( lb.x > lb.y ) 
  293.         l0 = lb.x;
  294.     else 
  295.         l0 = lb.y;
  296.     int rb;
  297.     if( l0 * 0.75 * M_SQRT2 + 1.0 == 1.0 ) 
  298.         rb = 0;
  299.     else
  300.         rb = (int)ceil(log4( M_SQRT2 * 6.0 / 8.0 * INV_EPS * l0 ) );
  301.     RecursivelyIntersect( a, 0., 1., ra, b, 0., 1., rb, parameters, index );
  302.     }
  303.     if( index < 9 )
  304.     parameters[0][index] = parameters[1][index] = -1.0;
  305.     Sort( parameters[0], index );
  306.     Sort( parameters[1], index );
  307.     return parameters;
  308.     }
  309.  
  310. void
  311. Bezier::ParameterSplitLeft( double t, Bezier &left )
  312.     {
  313.     left.p0 = p0;
  314.     left.p1 = new point( *p0 + t * ( *p1 - *p0 ) );
  315.     left.p2 = new point( *p1 + t * ( *p2 - *p1 ) ); // temporary holding spot
  316.     p2->refcount--;
  317.     p2 = new point( *p2 + t * ( *p3 - *p2 ) );
  318.     p1->refcount--;
  319.     p1 = new point( *(left.p2) + t * ( *p2 - *(left.p2) ) );
  320.     *(left.p2) = ( *(left.p1) + t * ( *(left.p2) - *(left.p1) ) );
  321.     (left.p3 = p0 = new point(*(left.p2) + t * (*(p1)-*(left.p2))))->refcount++;
  322.     left.p0->refcount++; left.p1->refcount++; 
  323.     left.p2->refcount++; left.p3->refcount++;
  324.     }
  325.     
  326. /*
  327.  * Intersect two curves, returning an array of two arrays of curves.
  328.  * The first array of curves corresponds to `this' curve, the second 
  329.  * corresponds to curve B, passed in.
  330.  * The intersection parameter values are computed by FindIntersections,
  331.  * and they come back in the range 0..1, using the original parameterization.
  332.  * Once one segment has been removed, ie the curve is split at splitT, the
  333.  * parameterization of the second half is from 0..1, so the parameter for
  334.  * the next split point, if any, must be adjusted.
  335.  * If we split at t[i], the split point at t[i+1] is 
  336.  * ( t[i+1] - t[i] ) / ( t - t[i] ) of the way to the end from the new
  337.  * start point.
  338.  */
  339.  
  340. Bezier **Bezier::Intersect( Bezier B )
  341.     {
  342.     // The return from FindIntersections will decrement all refcounts.
  343.     // (a c++-ism)
  344.     B.p0->refcount++; B.p1->refcount++; B.p2->refcount++; B.p3->refcount++; 
  345.     Bezier **rvalue = new Bezier *[2];
  346.     rvalue[0] = new Bezier[10];
  347.     rvalue[1] = new Bezier[10];
  348.     double **t = FindIntersections( *this, B );
  349.     int index = 0;
  350.     if( t[0][0] > -0.5 )
  351.     {
  352.     ParameterSplitLeft( t[0][0], rvalue[0][0] );
  353.     B.ParameterSplitLeft( t[1][0], rvalue[1][0] );
  354.     index++;
  355.     while( t[0][index] > -0.5 && index < 9 )
  356.         {
  357.         double splitT = (t[0][index] - t[0][index-1])/(1.0 - t[0][index-1]);
  358.         ParameterSplitLeft( splitT, rvalue[0][index] );
  359.         splitT = (t[1][index] - t[1][index-1])/(1.0 - t[1][index-1]);
  360.         B.ParameterSplitLeft( splitT, rvalue[1][index] );
  361.         index++;
  362.         }
  363.     }
  364.     rvalue[0][index] = *this;
  365.     rvalue[1][index] = B;
  366.     return rvalue;
  367.     }
  368.